Predicting criticality and dynamic range in complex networks: effects of topology 
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The collective dynamics of a network of coupled excitable systems in response to an external stimulus de- 
pends on the topology of the connections in the network. Here we develop a general theoretical approach to 
study the effects of network topology on dynamic range, which quantifies the range of stimulus intensities result- 
ing in distinguishable network responses. We find that the largest eigenvalue of the weighted network adjacency 
matrix governs the network dynamic range. Specifically, a largest eigenvalue equal to one corresponds to a crit- 
ical regime with maximum dynamic range. We gain deeper insight on the effects of network topology using a 
nonlinear analysis in terms of additional spectral properties of the adjacency matrix. We find that homogeneous 
networks can reach a higher dynamic range than those with heterogeneous topology. Our analysis, confirmed by 
numerical simulations, generalizes previous studies in terms of the largest eigenvalue of the adjacency matrix. 
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Numerous natural 01 01 and social 13] systems are accu- 
rately described as networks of interacting excitable nodes. 
The collective dynamics of such excitable networks often defy 
naive expectations based on the dynamics of the single nodes 
which comprise the network. For example, the collective re- 
sponse of a neural network can encode sensory stimuli which 
span more than 10 orders of magnitude in intensity, while 
the response of a single neuron (node) typically encodes a 
much smaller range of stimulus intensities. More generally, 
the range of stimuli over which a network's response varies 
significantly is quantified by dynamic range and is a funda- 
mental property, whether the network is comprised of people, 
cell phones^genes, or neurons. In neural networks, recent ex- 
periments 14j] suggest that dynamic range is maximized in a 
critical regime in which neuronal avalanches J5t] occur, con- 
firming earlier theoretical predictions J2]. It has been argued 
that this critical regime occurs when the effective mean 
degree of the network is one, i.e. the expected number of ex- 
cited nodes produced by one excited node is one. However, 
this criterion is invalid for networks with broad degree distri- 
butions A general understanding of how dynamic range 
and criticality depend on network structure remains lacking. 
In this Letter, we present a unified theoretical treatment of 
stimulus-response relationships in excitable networks, which 
holds for diverse networks including those with random, scale 
free, degree-correlated, and assortative topologies. 

As a tractable model of an excitable network, here we con- 
sider the Kinouchi-Copelli model J2], which consists of N 
coupled excitable nodes. Each node i can be in one of m 
states Xi. The state Xi = is the resting state, xi = 1 is 
the excited state, and there may be additional refractory states 
1. At discrete times t = 0, 1, ... the states 
of the nodes x\ are updated as follows: (i) If node i is in the 
resting state, x\ = 0, it can be excited by another excited node 
j, x l j = 1, with probability Ay, or independently by an ex- 
ternal process with probability rj. The network topology and 
strength of interactions between the nodes is described by the 
connectivity matrix A = {Aij}. In this model, rj is considered 



the stimulus strength, (ii) The nodes that are excited or in a re- 
fractory state, x\ > 1, will deterministically make a transition 
to the next refractory state if one is available, or otherwise re- 
turn to the resting state (i.e. x\ +1 = x\ + 1 if 1 < x\ < m — 1, 
anda^ +1 = Oif x\ =m- 1). 

An important property of excitable networks is the dynamic 
range, which is defined as the range of stimuli that is distin- 
guishable based on the system's response F. Following J2|, 
we quantify the network response with the average activity 
F = (f)t where (•)< denotes an average over time and /' is 
the fraction of excited nodes at time t. To calculate a system's 
dynamic range, we first determine a lower stimulus threshold 
rji ow below which the change in the response is negligible, and 
an upper stimulus threshold rjhigh above which the response 
saturates. Dynamic range (A), measured in decibels, is de- 
fined as A = 10 log 10 rjhigh/viow- To analyze the dynamics 
of this system, we denote the probability that a given node i is 
excited at time t by p\ . For simplicity, we will consider from 
now on only two states, resting and excited (m=2) jjj. Then, 
the update equation for p\ is 
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which can be obtained by noting that 1 — p\ is the probability 
that node i is resting at time t, and the term in large paren- 
theses is the probability that it makes a transition to the ex- 
cited state. We note that, in writing this probability, we treat 
the events of neighbors of node i being excited at time t as 
statistically independent. As noted before JHI^l!], this ap- 
proximation yields good results even when the network has a 
non-negligible amount of short loops. 

InRef. Jl, the response F was theoretically analyzed as 
a function of the external stimulation probability rj using a 
mean-field approximation in which connection strengths were 
considered uniform, Aij = a/N for all It was shown that 
at the critical value a = 1, the network response F changes 
its qualitative behavior. In particular, lim F — if a < 1 
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and lim F > if a > 1. In addition, the dynamic range of 

the network was found to be maximized at a = 1. The pa- 
rameter a is defined in Refs. J2, 01 as an average branching 
ratio, written here as a = ^X^j^'j = (d m ) = (d° ut ), 
where d" 1 ^ A ■ ■ ""^ ^ otl * 



X , A„ and d° 



J2j Aji are the in- and 



out-degrees of node i, respectively, and (•} is an average over 
nodes. For the network topology studied by Ref. [0] cr = 1 
marks the critical regime in which the expected number of 
excited nodes is equal in consecutive timesteps. Such criti- 
cal branching processes result in avalanches of excitation with 
power-law distributed sizes. Cascades of neural activity with 
such power-law size distributions have been observed in brain 
tissue cultures J3l, awake monkeys (12], and anesthetized rats 
1 1311 . While cr=l successfully predicts the critical regime 
for Erdos-Renyi random networks |0], this prediction fails 
in networks with a more heterogeneous degree distribution 
10,0]. Perhaps more importantly, previous theoretical anal- 
yses 0, do not account for features that are commonly 
found in real networks, such as community structure, correla- 
tion between in- and out-degree of a given node, or correlation 
between the degree of two nodes at the ends of a given edge 
1 1411 . Here, we will generalize the mean-field criterion a = 1 
to account for complex network topologies. 

To begin, we note that lim F = corresponds to the 
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fixed point p = of Eq. (Q]i with 77 = 0. To examine the 
linear stability of this fixed point, we set 77 = and lin- 
earize around = 0, assuming p\ to be small, obtaining 
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Thus, the stability of the solution p = is governed by the 
largest eigenvalue of the network adjacency matrix, A, with 
A < 1 being stable and A > 1 being unstable. Therefore, 
the critical state described in previous literature, occurring at 
various values of (d), should universally occur at A = 1. Im- 
portantly, since Aij > 0, the Perron-Frobenius theorem guar- 
antees that A is real and positive 1 1 511 - Other previous studies 
in random networks have also investigated spectral properties 
of A to gain insight on the stability of dynamics in neural net- 
works 1 161] and have shown how A could be changed by modi- 
fying the distribution of synapse strengths lfl7ll . An important 
implication of Eq. (O is that, when p and rj are small enough, 
p should be almost proportional to the right eigenvector u cor- 
responding to A, so we write pi — Cui + ei, where C is a 
proportionality constant and the ej error term captures the de- 
viation of actual system behavior from the linear analysis. To 
first order, the constant C is related to the network response F 
since, neglecting e, we have 

i i 

The linear analysis allowed us to identify A = 1 as the point 
at which the network response becomes non-zero as 77 — >• 0. 
In what follows, we use a weakly nonlinear analysis to obtain 
approximations to the response F(rj) when 77 is small. As we 



will show, these approximations depend only on a few spec- 
tral properties of A. Assuming AijPj <C 1 (which is valid 
near the critical regime if each node has many incoming con- 
nections), we approximate the product term of Eq. (fl} with an 
exponential, obtaining in steady state 



Pi = (1 - Pi) [v + (1 - v) 
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which we expand to second order using Eq. (O and Ai 
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To eliminate the error term ej from Eq. (0, we multiply by Vi, 
the ith entry of the left eigenvector corresponding to A, and 
sum over i. We use the fact that v T Ae = Xv T e, where v T de- 
notes the transpose of v, and neglect the resulting small term 
(1 — A) J2i v i e i close to the critical value A = 1, obtaining 



C(uv) = r)((v)-C(uv))+(l-r))CX(uv)- 
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This equation is quadratic in C [and therefore in F, via 
Eq. and linear in rj, and may be easily solved for either. 
For -q = the nonzero solution for F is 

p - (A-i) HW m 

*"=°- (A + |A 2 ) (u*v) ■ (/) 

A more refined approximation than Eq. © can be obtained by 
repeating this process without expanding Eq. (|4J), which yields 
the linear equation for 77 

C{uv) =53(l-Cti i )(»7 + (l-77)[l-exp(-AC« i )]). (8) 
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Before numerically testing our theory, we will explain how 
it relates to previous results. For a network with correlations 
between degrees at the ends of a randomly chosen edge (assor- 
tative mixing by degree lfl4ll ). measured by the correlation co- 
efficient p = (d™d° ut ) e /(d in d otl *), with (-} e denoting an av- 
erage over edges, the largest eigenvalue may be approximated 
by A w p(d m d out )/(d) Q. In the absence of assortativity, 
when p = 1, A w (d m d out ) / (d) . If, in addition, there are no 
correlations between di„ and d out (node degree correlations) 
or if the degree distribution is sufficiently homogenous, then 
(d m d out ) ps (d) 2 and the approximation reduces to A ps (d). 
In the case of Ref. 0], A ps (d) applies, and in the case of 
Refs. |@,0, A ps (d m d out )/(d) applies. 

We test our theoretical results via direct simulation of the 
Kinouchi-Copelli model on six categories of directed net- 
works with N = 10, 000 nodes: (category 1) Random net- 
works with no node degree correlation between d m and d out ; 
(category 2) Random networks with maximal degree correla- 
tion, d m = d out ; (category 3) Random networks with moder- 
ate correlation between d m and d out ; (category 4) Networks 
with power law degree distribution with power law exponents 
7 G [2.0,6.0], with and without node degree correlations; 
(category 5) Networks constructed with (d) = 1, and assorta- 
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tivity coefficient p varying in [0.7, 1.3]; (category 6) Networks 
with weights which depend on the degree of the node from 
which the edge originates, A^ = a/d° ut . 

We created networks in multiple steps: first, we created 
binary networks {A^j e {0, 1}) with target degree distribu- 
tions as described below; next, we assigned a weight to each 
link, drawn from a uniform distribution between and 1 ; fi- 
nally, we calculated A for the resulting network and multi- 
plied A by a constant to rescale the largest eigenvalue to the 
targeted eigenvalue. This process was restarted from the first 
step for every network used in categories 1-4, creating a struc- 
turally different network for each simulation. The initial bi- 
nary networks in categories 1-3 were Erdos-Renyi random 
networks, constructed by linking any pair of nodes with prob- 
ability p = 10 /N iflill . Maximal degree correlation resulted 
from creating undirected binary networks and then forcing 
Ay = Aji for i < j while assigning weights. Moderate de- 
gree correlation resulted from making undirected binary net- 
works but allowing Aij ^ Aji when weights were assigned. 
The algorithms for constructing the initial binary networks 
of categories 4-6 placed links randomly between nodes with 
specified in- and out-degrees via the configuration model lfl9ll . 
For this model, we generated in- and out-degree sequences 
from a power law distribution of desired exponent 7 by calcu- 
lating the expected integer number of nodes with each integer 
degree, from minimum degree 10 to maximum degree 200. In 
creating category 5 networks, we initially created one scale 
free network with power law exponent 7 = 2.5 and A = 1. 
Then, to change the degree of assortativity, we modified this 
original network by choosing two links at random and swap- 
ping them if the resulting swap would change the assortativity 
in the direction desired. This process was repeated until a 
desired value of p was achieved. Importantly, this swapping 
makes it possible to leave the degree distributions of the net- 
work unchanged, while still changing the assortative or disas- 
sortative properties of the network as in lfl4l Eoll . Therefore, 
by this method we may maintain exactly the same degree dis- 
tribution and mean degree, yet modify A by virtue of A oc p. 

In the six network types tested, results of simulations unan- 
imously confirm the hypothesis that criticality occurs only for 
largest eigenvalue A = 1. We present representative results 
in Fig. Q] (a), noting that each line and set of points corre- 
sponds to a single network realization, implying that the effect 
of the largest eigenvalue on criticality is robust for individual 
systems. Fig. Q] (a) shows the response F as a function of 
stimulus r\ for scale-free networks with exponent 7 = 2.5, 
constructed with no correlation between in- and out-degree, 
highlighting the significant difference between the regimes of 
A < 1 and A > 1, with the critical data corresponding to 
A = 1. The lines were obtained by using Eqs. d3) and ((8). 
Fig. [T](b) shows A as a function of A, using r\high = 1 an d 
Viow = 0.01, with the maximum occurring at A = 1. Sim- 
ilar results showing criticality and maximum dynamic range 
at A = 1 are obtained for networks of all categories 1-5. Fig. 
|2]shows F v ^q for networks of categories 3-5, confirming the 
transition predicted by the leading order analysis in Eq. (|2j. 
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FIG. 1 : (color online) (a) Response F vs. stimulus r\ for power law 
networks with exponent 7 = 2.5 and no correlation between di n and 
dout. Eq. (lines) captures much of the behavior of the simulation 
(circles), particularly for low levels of r) and F, as expected from 
approximating Eq. (TJ. (b) Dynamic range A is maximized at A = 1 
in both simulation results (circles) and Eq. l|6j (line). 

The symbols show the result of direct numerical simulation 
of the Kinouchi-Copelli model, the solid lines were obtained 
by iterating Eq. (OJ, and the dashed lines were obtained from 
Eq. ©■ Fig. 12 a) shows that criticality occurs at A = 1 (indi- 
cated by a vertical arrow) rather than at (d) = 1 for a category 
3 random network. Fig. [2b) shows that criticality occurs at 
A = 1 for scale-free networks (category 4). Correlations be- 
tween di n and d out affect the point at which A = 1 occurs 
(vertical arrows). In Fig. |2 C X the mean degree was fixed at 
(d) = 1, while A was changed by modifying the assortative 
coefficient p. As predicted by the theory, there is a transition 
at A = 1 even though the mean degree is fixed. 

We now explore the question of what network topology will 
best enhance dynamic range. In many of the systems we sim- 
ulate, a majority of the variation in dynamic range from one 
stimulus-response curve to another occurs due to variation 
at the low stimulus end of the curve, since most of the sys- 
tems tend to saturate at around the same high stimulus levels 
(though this may not be the case for neuronal network exper- 
iments |0]). We therefore consider the following approximate 
measure of dynamic range, A, obtained by setting r\high to 
one in the definition of A, A = 10 log 10 1 /r)*, where 77* is the 
stimulus value corresponding to a lower threshold response 
F*. Since dynamic range is maximized at criticality, we set 
A = 1, solve Eq. (O for substitute it into the definition of 
A using Eq. (O, retaining the leading order behavior to get 

2 (vu 2 ) 
Km ax = 10 log 10 - 10 log 10 j^j-fe. (9) 
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(a) random, degree correlated 



(b) scale free 



(c) scale free, assortative 
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FIG. 2: (color online) i^_).o obtained from direct numerical simulation of the Kinouchi-Copelli model (symbols) plotted against (d) (a, b) 
and A (c). Blue solid lines result from iterating Eq. {T} and green dashed lines result from Eq. Q. Small arrows show where A = 1 predicts a 
phase transition, (a) A set of random networks (category 3) showing that criticality occurs at A = 1 (arrow), but not (d) = 1. (b) Criticality 
in scale free networks (category 4) with node degree correlation also occurs at A = 1 (arrow), but not (d) = 1. (c) Category 5 networks are 
tuned through criticality by changing assortativity, without changing the degree distributions and fixed (d) = 1. 



The first term of this equation shows that Am ax depends 
on F*. Since the entries of the right (left) dominant eigen- 
vector are a first order approximation to the in-degree (out- 
degree) of the corresponding nodes 12111 . the second term 
suggests that maximum dynamic range should increase (de- 
crease) as the degree distribution becomes more homogenous 
(heterogeneous). For example, consider the case of an undi- 
rected, uncorrelated network, in which t>, = it, ~ dj . The sec- 
ond term is then approximately —10 log 10 ((d 3 )/ (d) 3 ), which 
is maximized when di is independent of i. This corroborates 
the numerical findings in Refs. IH 13] that random graphs 
enhance dynamic range more than more heterogeneous scale 
free graphs, and that the heterogeneity of the degree distribu- 
tion affects dynamic range Q]. To test our result, we simu- 
late scale free networks with different power law exponents 
7 £ [2.0, 6.0], yet with A = 1 to maximize dynamic range in 
each case. Results of simulation (circles) plotted against the 
prediction of Eq. (O (line) are shown in Fig. [3] 

In summary, we analytically predict and numerically con- 
firm that criticality and peak dynamic range occur in networks 
with largest eigenvalue A = 1. This result holds for diverse 
network topologies including random, scale-free, assortative, 
and/or degree-correlated networks, and for networks in which 
edge weights are related to nodal degree, thus generalizing 
previous work. Moreover, we find that homogeneous (het- 
erogeneous) network topologies result in higher (lower) dy- 
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FIG. 3: For power-law degree distributions with A = 1, peak dy- 
namic range increases monotonically with network homogeneity, as 
measured by power law exponent 7. Simulations (circles) agree well 
with our predictions [Eq. line]. 



namic range. Previous demonstrations of how A governs net- 
work dynamics in many other models (see 12111 and references 
therein) suggest that the generality of our findings may extend 
beyond the particular model studied here. Previous model 
studies have shown that mutual information between stimulus 
and response is also maximized at criticality @]. Our findings 
suggest that peak mutual information will also be determined 
by A = 1, but verifying this will require additional investi- 
gation. Taken together with related experimental findings J3l, 
our results are consistent with the hypotheses that 1) real brain 
networks operate with A « 1, and 2) if an organism benefits 
from large dynamic range, then evolutionary pressures may 
act to homogenize the network topology of the brain. 
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